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ABSTRACT 


This  study  explores  the  fluxes  of  heat  and  salt 
associated  with  thermohaline  staircases  in  the  Beaufort  Sea. 
An  inverse  model  is  developed  to  calculate  the  vertical 
transport  of  properties  in  the  southern  portion  of  the 
Beaufort  Sea  directly  from  observations.  The  applicability 
of  laboratory  derived  4/3  flux  law  is  addressed.  Three 
formulations  of  the  static  advective-dif fusive  equations  are 
discretized  on  a  uniform  grid,  and  inverted  using  the  method 
of  Total  Least  Squares.  The  first  formulation  is  based  on 
the  classic  flux-gradient  model,  and  is  analyzed  in  both  one 
and  three-dimensions.  The  second  formulation  utilizes 
Turner  formulation  of  the  4/3  flux  law,  but  the  formulation 
of  the  coefficient  C  (Rp)  remains  unknown.  The  third 
formulation  includes  the  full  form  of  the  Kelley  model,  but 
the  amplitude  of  C (Rp)  is  unknown.  Layer  averaged  heat  flux 
through  the  thermohaline  staircases  is  found  to  be  on  the 
order  of  1  W/m^,  which  is  significantly  higher  than  previous 
studies.  This  implies  that  the  laboratory  formulation  of  the 
4/3  flux  law  may  require  calibration  in  order  to  accurately 
represent  Arctic  conditions. 
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I.  INTRODUCTION  MUD  BACKGROUND 

A.  ARCTIC  CIRCULATION 

Sea  ice  in  the  Arctic  is  melting  at  an  alarming  rate, 
much  faster  than  expected  (Stroeve  et  al . ,  2007) .  The  years 
2007-2009  have  seen  the  lowest  summer  sea  ice  coverage  in 
the  Arctic  in  the  modern  era.  Ice  melt  in  the  Arctic  has 
been  largely  attributed  to  rising  greenhouse  gas 
concentrations  and  the  associated  warming  of  the  atmosphere. 
The  ocean,  however,  responds  to  the  atmospheric  conditions 
on  much  longer  time  scales.  We  are  only  beginning  to 
understand  the  response  of  the  Arctic  Ocean  to  the 
atmospheric  forcing  associated  with  climate  change  (Stroeve 
&  Maslowski,  2008;  Dickson  et  al . ,  2000) .  One  such  response 
is  the  delivery  of  significantly  elevated  quantities  of 
warmer  water  to  the  Canada  Basin  and  Beaufort  Sea.  This 
additional  warm  water,  which  has  taken  decades  to  traverse 
the  Canada  Basin,  adds  significant  amounts  of  heat  to  the 
upper  Arctic  Ocean,  aggravating  the  ice  melt  problem. 

The  heat  stored  in  the  oceans  directly  influences  the 
melting  rate  of  the  ice  from  below.  Perovich  et  al .  (2003) 

showed  a  1-2  W/m^  difference  between  the  amount  of  heat 
necessary  to  melt  the  ice  from  below  and  the  atmospheric 
heat  fluxes  calculated  over  open  water.  They  suggest  that 
there  must  be  oceanic  sources  of  heat  to  the  ocean-ice 
interface.  Understanding  the  mechanisms  that  drive  the 

heat  flux  to  the  ice  from  the  ocean  is  critical  to  the 
explanation  and  prediction  of  increased  sea  ice  melt.  The 
ocean  can  play  a  role  in  melting  sea  ice  via  several 
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mechanisms.  Ice-ocean  interactions  include  horizontal 
advection  of  warm  water  from  the  Pacific  and  Atlantic  oceans 
under  the  ice  cove;  the  accumulation  of  heat  due  to 
increased  solar  radiation  in  ice  diminished  regions;  and 
locally  induced  upward  heat  flux  into  the  mixed  layer  due  to 
upwelling,  topographically  controlled  flow,  and  eddies 
(Maslowski  et  al .  ,  2009). 

This  thesis  will  focus  on  the  role  of  double  diffusion 
in  oceanic  heat  transfer,  which  is  a  subject  of  growing 
interest  in  oceanography  and  Arctic  sciences.  It  is  argued 
that  double-diffusive  convection  may  also  play  a  significant 
role  in  transferring  heat  upward  into  the  Arctic  mixed 
layer.  The  heat  transported  toward  the  ocean  surface  via 
diffusive  convection  may  be  a  critical  component  to 
understanding  the  connection  between  ocean  dynamics  and  sea 
ice  extent  in  the  Arctic. 

B.  DOUBLE-DIFFUSIVE  CONVECTION 

Double  diffusion  is  defined  as  a  set  of  processes 
related  to  the  difference  in  molecular  dif fusivities  of  heat 
and  salt.  It  is  known  to  influence  vertical  mixing  in  the 
ocean,  and  comes  in  two  forms:  salt  fingering  and  diffusive 
convection.  Salt  fingering  is  prevalent  in  the  mid¬ 
latitudes,  while  diffusive  convection  is  the  form  of  double 
diffusion  seen  in  the  Arctic. 

The  term  diffusive  convection  is  used  to  describe  two 
types  of  motion  that  arise  in  fluids  with  a  negative 
vertical  gradient  of  temperature  and  salinity.  The  first  is 
characterized  by  oscillations  in  smooth  temperature-salinity 
gradients.  The  other,  diffusive  layering,  operates  in  a 
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thin  interface,  which  separates  two  homogeneous  layers.  Heat 
and  salt  are  transferred  across  the  diffusive  interface. 


The  oscillatory  regime  of  double  diffusive  convection 
is  described  in  Figure  1 . 


cold,  fresh 


heat  ‘ '  heat 

1=^  1  <=1 


warm,  salty 


Figure  1.  The  oscillatory  regime  of  double  diffusive 

convection 


Imagine  a  parcel  of  water  in  a  lower  layer  that  is 
displaced  upward.  This  parcel  is  warmer  and  more  saline 
than  the  layer  of  water  above.  As  the  parcel  enters  the 
layer  above,  it  immediately  begins  to  lose  heat,  via 
molecular  diffusion,  to  the  surrounding  upper  layer  at  a 
rate  nearly  100  times  as  fast  as  it  loses  salt.  The  result 
is  cooling  of  the  parcel,  while  it  maintains  its  salinity. 
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The  parcel  then  begins  to  descend,  due  to  its  increasing 
density,  as  it  cools.  The  denser  parcel  begins  to  gain  heat 
in  the  lower,  warmer  layer,  again  via  molecular  diffusion, 
and  once  it  gains  a  sufficient  amount  of  heat  its  density 
decreases  and  the  parcel  begins  to  rise  again.  This  mixing 
results  in  the  formation  of  layers  of  near  constant 
temperature  and  salinity.  It  is  these  layers  that  define 
the  second  form  of  diffusive  convection,  diffusive  layering. 
This  process  is  known  to  be  prevalent  in  the  Beaufort  Sea. 


Potential  Temperature  ( C-) 


Figure  2.  Temperature  -  Salinity  (From  Timmermans,  2008) 


Diffusive  layering  occurs  in  a 
where,  in  the  case  of  the  Arctic,  cold 
warm  and  salty  water,  as  seen  in  Figure 


stable  environment, 
fresh  water  overlies 

2  . 
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Figure  3.  Diffusive  Layering  -  Fluxes  of  heat  and  salt  are 
transported  upward  across  the  thin  layer,  causing 
convection  and  mixing  within  the  homogeneous  layers. 


Diffusive  layering  is  manifested  in  the  Arctic 
thermocline  by  a  series  of  homogeneous  layers  of  near 
constant  temperature  and  salinity  separated  by  a  thin 
diffusive  interface.  This  natural  occurrence  is  often 
referred  to  as  a  thermohaline  staircase.  The  physical 
processes  associated  with  the  maintenance  of  this  structure 
(Figure  3)  can  be  described  as  follows:  Heat  and  salt  are 
transferred  upward  through  the  diffusive  interface.  The 
water  just  above  the  interface  gains  sufficient  heat  to 
begin  upward  motion  due  to  excess  of  buoyancy.  This  water 
begins  to  rise  within  the  layer,  driving  the  convection  that 
ultimately  maintains  the  homogeneous  properties  of  the  mixed 
layer . 
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The  data  shown  in  Figure  2  are  representative  of  the 
majority  of  the  data  collected  by  the  Woods  Hole 
Oceanographic  Institute  (WHOI)  Ice  Tethered  Profiler  (ITP) 
program.  The  warm  water  seen  from  ~300  meters  and  below 


originates  in 

the 

Pacific 

and 

Atlantic 

oceans 

and  is 

advected  into 

the 

Beaufort 

Sea 

via  the 

Arctic 

current 

system.  The 

cold 

water  above 

is  a  result  of 

downward 

temperature  flux  associated  with  the  cold  surface  and  air- 
sea-ice  interaction.  The  less  saline,  or  relatively  fresh, 
surface  water  is  a  result  of  annual  ice  melt  and  fresh  water 
flux  associated  with  river  runoff.  These  conditions,  and 
the  lack  of  vertical  mixing  below  the  surface  mixed  layer, 
provide  an  ideal  environment  for  the  formation  of  double 
diffusive  staircases. 


C.  EARLIER  ESTIMATES  OF  HEAT  FLUX 


The  fluxes  associated  with  double  diffusion  have  been 
studied  extensively.  It  is  generally  accepted  that  in 
systems  characterized  by  active  diffusive  layering,  the 
interaction  between  layers  is  controlled  by  the  flux  of  salt 
and  heat  across  the  diffusive  interfaces.  Based  on  a  set 
laboratory  experiments  the  balance  of  temperature  variation 
between  adjacent  layers  and  heat  flux  (Turner,  1973)  was 
proposed  as  follows. 
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aFj.  =  Ay-i^aATy  (1) 


aAT 


(2) 


This  relationship  (1)  has  since  been  known  as  the  4/3 
flux  law.  Ai  has  the  dimensions  of  velocity  and  is  a 
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function  of  the  density  ratio  (Rp)  .  Ft  is  temperature  flux, 
a  is  the  thermal  expansion  coefficient,  and  (3  is  the  saline 
contraction  coefficient. 

The  4/3  flux  law  states  that  the  fluxes  associated  with 
double  diffusion  are  controlled  entirely  by  the  variation  in 
temperature  or  salinity  across  the  diffusive  interface  and 
is  independent  of  layer  thickness.  The  lack  of  sensitivity 
of  fluxes  to  vertical  scales  is  critically  important  to  our 
observation-based  study  because  the  thickness  of  an 
interface  is  often  much  smaller  (<25cm)  than  the  resolution 
of  the  sensor,  whereas,  a  difference  in  mean  temperatures  of 
each  layer  can  easily  be  measured.  Dimensional  analysis 

suggests  that  the  4/3  flux  law  can  be  written  as  follows: 

1 

aF^=C[R^)\^^  ■{a^Tf  (3) 

V  ^  / 

C  (Rp)  is  an  unknown  constant,  g  =  9.8  m/s,  v  =  1.8X10”® 
m^/s  is  the  kinematic  viscosity,  k  =  1.4X10”^  m^/s  is  the 

molecular  diffusivity  of  heat. 

Marmorino  and  Caldwell  (1976)  examined  the  4/3  flux  law 
in  depth.  They  discovered  an  empirical  formulation  of  the 
coefficient  C (Rp)  using  the  following  equations: 

=  (4) 

-^  =  0.101exp|4.6exp[^-0.54(^i?^ -l)  |  (6) 

^  sp 
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c(i?^)  =  0.00859 exp (4.6 exp [-0.54 -l}]) 


(7) 


Using  Huppert's  formula  (4),  they  were  able  to 
determine  a  relationship  between  the  measured  heat  flux  (H) 
and  (5) ,  the  theoretical  heat  flow  through  a  non-def ormable 
interface  (Hgp)  .  They  proposed  a  new  formulation  also  based 
on  laboratory  experiments  (6),  which  ultimately  yielded  a 
value  for  C  (Rp)  . 


The  coefficient  C (Rp)  has  since  been  scrutinized,  most 
notably  by  Kelley  (1990)  .  In  his  1990  paper,  Kelley 
proposed  a  new  value  for  C  (Rp)  based  on  what  he  referred  to 
as  a  full  collection  of  laboratory  measurements.  He 
analyzed  laboratory  experiments  for  l<i?^<10,  while  Turner 
(1973)  only  looked  at  values  of  Rp  up  to  7. 


C(i?^)  =  0.0032  exp 


4.8 


0.72 


y 


(8) 


Equation  (8)  is  the  result  of  Kelley's  laboratory 
experiments,  and  along  with  (7),  the  two  have  been  used 
interchangeably  as  a  standard  method  for  calculating  fluxes 
using  the  4/3  flux  law.  However,  Kelley  suggested  that 
perhaps  the  4/3  exponent  is  incorrect  and  should  be  changed 
to  5/4.  Additional  experimentation  is  required  to  support 
this  proposition. 

Padman  and  Dillon  (1987)  showed  that  using  the 
Marmorino  and  Caldwell  formulation  of  the  4/3  flux  law,  they 
could  calculate  the  fluxes  in  the  Beaufort  Sea.  Their 
results  showed  that  those  fluxes  were  0.02  —  0.1  Wm”^ .  A 
more  recent  study  (Timmermans,  2008),  also  based  on 
extrapolation  of  laboratory  derived  results,  estimated  the 
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heat  flux  associated  with  the  diffusive  staircases  as 
measured  by  the  Ice  Tethered  Profilers  to  be  0.22  +/-  0.10 
W/m^ . 

Wilson  (2007)  compared  the  two  formulations  (7)  and  (8) 
using  data  collected  via  the  Ice  Tethered  Profilers,  a  set 
of  moored  CTDs  in  the  Beaufort  Sea,  which  will  be  discussed 
in  detail  in  Chapter  II.  She  found  that  in  order  to  fit  the 
data  collected  in  the  Beaufort  Sea  to  these  formulae  each 
needed  to  be  multiplied  by  a  transfer  coefficient.  She 
suggested  that  the  laboratory  results  of  both  Kelley  (1990) 
and  Marmorino  and  Caldwell  (1976)  cannot  be  applied  to 
actual  data  collected  in  the  field,  and  perhaps  the 
laboratory  measurements  require  calibration  prior  to 
application  in  the  Arctic.  She  concluded  that  the  average 
heat  flux  in  the  Beaufort  Gyre  in  excess  of  1  Wm”^ .  This  is 
significantly  larger  than  estimates  using  the  original  4/3 
flux  law. 

Caro  (2009)  used  a  set  of  high-resolution  numerical 
experiments  to  validate  some  of  the  conclusions  made  by 
Wilson.  He  used  two-dimensional  numerical  models  to 
estimate  the  heat  fluxes  in  the  Beaufort.  His  efforts 
resulted  in  fluxes  on  the  order  of  1  W/m^. 

Timmermans'  findings  differ  from  Wilson's  (2007)  and 
Caro's  (2009)  by  an  order  of  magnitude,  and  in  some  cases 
Padman  and  Dillon's  differ  by  two  orders  of  magnitude.  Such 
a  wide  range  of  estimates  in  a  region  critically  important 
for  climate  variability  provides  the  impetus  for  this  study. 
This  thesis  attempts  to  resolve  the  controversy  in 
assessment  of  the  vertical  heat  transport  by  applying 
inverse  modeling  techniques  to  observations. 
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D .  INVERSE  MODELING 

The  study  of  the  oceans  traditionally  follows  two 
distinct  directions.  One  involves  numerical  modeling  and 
theoretical  manipulation  of  the  equations  governing  dynamics 
and  thermodynamics  of  the  sea.  The  other  is  observationally 
oriented,  where  real  measurements  are  discussed  in  the 
context  of  qualitative  physical  principles. 

As  oceanography  matures  into  a  quantitative  science,  it 
becomes  increasingly  desirable  to  connect  the  observational 
and  modeling  aspects  of  the  field.  Laboratory  experiments 
often  attempt  to  bridge  this  gap  between  hard  facts  and 
theoretical  deductions.  Laboratory  researchers  strive  to 
replicate  the  ocean  environment  as  closely  as  possible,  but 
exact  correspondence  is  seldom  achieved,  given  the 
substantial  differences  in  the  conditions  and  spatial  scales 
between  the  laboratory  and  nature. 

Another  more  recent  approach,  which  makes  it  possible 
to  connect  data  with  theory,  involves  inverse  modeling. 
Inverse  models  attempt  to  directly  calculate  model 
parameters  from  observations  by  minimizing  errors  associated 
with  observational  measurements  and  fitting  those 
measurements  to  the  governing  equations.  In  our  case, 
temperature,  salinity  and  depth  measurements  were  taken  by 
the  Ice  Tethered  Profilers,  and  velocity  and  diffusion 
coefficients  are  calculated  from  those  observations  using 
the  method  of  Total  Least  Squares. 
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(11) 


V*v  +  — =  0 

dz 

The  inverse  model  in  this  study  is  based  on  the  static 
advective  diffusive  (9),  (10)  and  mass  continuity  equations 

(11)  in  divergence  form.  Ft,  and  Fs  is  the  temperature  or 
salinity  flux.  The  formulation  of  Ft  and  Fs  is  central  to 

this  study.  In  this  work  we  examine  flux-gradient 

formulations  in  which  Ft,  and  Fs  are  assumed  to  be 
controlled  by  the  background  temperature  and  salinity 

gradients,  as  well  as  the  interfacial  flux  laws.  Results  of 
both  are  qualitatively  consistent,  but  differ  in  details. 


E. 


UNRESOLVED  ISSUES 


This  thesis  attempts  to  quantify  the  heat  flux 
associated  with  diffusive  convection  independent  of  the 
extrapolation  of  laboratory  derived  laws.  This  information 
is  used  to  address  the  following  questions. 

1,  What  are  the  typical  vertical  heat  and  salt  fluxes 

in  the  diffusive  staircases  in  the  Beaufort  Gyre?  Vertical 
heat  fluxes  of  0.02  —  6  Wm“^  represent  an  extremely  wide 
range  making  it  difficult  to  assess  the  large-scale 
consequences  of  heat  transfer  in  the  Arctic.  If  the  value 
is  less  than  1  Wm"^,  then  it  is  likely  that  the  heat 

delivered  upward  vis-a-vis  diffusive  convection  is  not  a 
major  contributor  to  sea  ice  melt.  However,  if  Wilson 
(2007)  is  correct,  then  diffusive  convection  is  sufficient 
to  affect  the  pattern  on  stratification  and  ocean  climate  in 
the  Arctic. 

2.  Can  the  heat  flux  through  the  diffusive  staircases 
significantly  contribute  to  the  heat  imbalance  associated 
with  the  melting  of  the  polar  ice  cap? 
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3.  What  are  the  relative  roles  of  the  vertical  and 
lateral  transport  of  heat  in  Beaufort  gyre?  Is  it  possible 
to  capture  the  zero-order  dynamics  at  play  by  a  one¬ 
dimensional  horizontally  averaged  model? 

4.  Is  it  possible  to  accurately  evaluate  fluxes 
associated  with  double  diffusion  by  inverse  modeling 
techniques?  Lee  and  Veronis  (1991)  showed  that  it  is 
possible  to  quantify  the  fluxes  via  an  inverse  model  using 
the  method  of  Total  Least  Squares  in  the  tropical  Atlantic 

Ocean.  a  similar  inverse  model  of  the  advective-dif fusive 
equations  is  applied  to  data  collected  from  2004-2009  via 
the  Ice  Tethered  Profilers  operated  by  the  Woods  Hole 
Oceanographic  Institute.  Using  methods  described  in  Lee  and 
Veronis  (1991)  the  inverse  technique  is  implemented  to 
determine  eddy  diffusivity  coefficients  using  the  method  of 
total  least  squares.  A  comparison  to  the  estimates  using 
the  4/3  flux  laws  described  in  Timmermans  (2008)  and  Padman 
and  Dillon  (1987)  is  discussed  along  with  a  comparison  to 
Wilson  (2007),  additionally  the  results  of  this  study  area 
compared  to  the  model  results  from  Caro  (2009) . 

This  thesis  is  organized  as  follows:  Chapter  II  reviews 
the  data  from  the  Ice  Tethered  Profilers  and  some  of  the 
unique  processing  required  for  obtaining  the  appropriate 
data  set  for  the  inverse  calculation.  Chapter  III  offers  a 
detailed  description  of  the  mathematical  foundation  of  the 
inverse  model  as  well  as  its  application  to  the  ITP  data 
set.  Chapter  IV  presents  the  results  of  four  models,  the 
one  dimensional  (vertical  only)  model,  which  uses  an  average 
single  profile  to  determine  fluxes,  and  three  variations  of 
a  full  three  dimensional  model.  A  discussion  of  the 
findings  and  conclusions  is  in  Chapter  V. 
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II.  DATA  AND  DATA  PROCESSING 


A.  ICE-TETHERED  PROFILER  DATA 

The  Woods  Hole  Oceanographic  Institute  (WHOI)  Ice 
Tethered  Profiler  (ITP)  program  was  established  in  2004  to 
provide  timely  and  accurate  oceanographic  measurements  in 
the  Arctic  Ocean.  Historically,  access  to  the  ice  covered 
Arctic  Ocean  has  been  limited.  The  ITPs  assist  in  assuring 
consistent  measuring  capability  in  a  data-sparse 
environment.  Between  2004  and  2006,  the  first  six  ITPs  were 
deployed  in  the  Beaufort  Sea.  Data  from  ITPs  1-6  are  used 
in  this  study.  As  of  July  2009,  these  six  sensors  collected 
over  7000  profiles,  all  of  which  have  been  made  available  to 
the  public. 

The  ITP  is  a  moored  conductivity,  temperature,  and 
depth  (CTD)  sensor.  Each  one  is  designed  to  move  with  the 
ice  it  is  moored  to  and  has  the  potential  lifespan  of 
approximately  three  years.  The  ITP  consists  of  a  small 
platform  on  top  of  which  rests  a  buoy.  A  10-inch  hole  is 
drilled  through  the  sea  ice  and  the  instrument  is  inserted 
in  the  water  from  above.  The  instrumentation  is  tethered  to 
the  surface  float  via  a  wire  rope  that  is  weighted  on  the 
bottom.  This  rope  is  approximately  800  meters  in  length, 
and  the  instrument  traverses  the  wire  rope  several  times 
daily . 
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Figure  4 
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Figure  5.  Location  of  ITP  1-6 


Figure  4  illustrates  the  ITP  system  and  Figure  5 
presents  the  trajectories  of  the  first  six  profilers 
followed  over  their  lifespans.  The  data  used  in  this  study 
are  labeled  "Level  3  Archive  Data"  by  WHOI .  This  form  of 
data  represents  the  best  possible  representation  of  the 
ocean  properties  as  measured  by  the  ITPs.  Corrections  have 
been  applied  that  account  for  sensor  response,  regional 
conductivity  variation,  and  quality  assurance.  Each  data 
file  is  formatted  to  be  compatible  with  MAT  LAB .  Each 
profile  was  retrieved  via 
f tp : / /f tp . who! . edu/whoinet / itpdata/  and  processed  further  as 
described  below. 
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B.  DATA  PROCESSING 

The  corrected  ITP  data  described  above  was  arranged  by 
profiler  number  1-6  and  processed  together  as  one  data  set. 
Processing  began  with  the  corrected  ITP  data  obtained  from 
the  WHOI  ITP  Web  site.  Wilson  (2007)  suggested  that 
anomalous  results  were  obtained  from  the  profiles 
corresponding  to  the  downward  data  collection  for  each  cast. 
This  is  possibly  due  to  sensor  contamination  associated  with 
the  main  ITP  sensor  platform  below  the  sensor  suite.  As  the 
profiler  is  lowered  into  the  water  the  platform  interacts 
with  the  water  just  below  the  sensors,  affecting  the 
measurements.  Based  on  this  suggestion,  all  downward 
collected  profiles  were  removed  from  the  total  dataset, 
effectively  reducing  the  amount  of  data  by  half. 

Diffusive  convection  occurs  in  the  thermocline  region 
of  the  Beaufort  Sea  between  approximately  200  meters  and  400 
meters  depth,  with  the  majority  of  the  diffusive  layers  in 
the  upper  100  meters  of  the  thermocline.  The  data  were  then 
reduced  to  only  those  data  from  200  db  pressure  to  300  db 
pressures,  which  correspond  to  roughly  195  meters  to  295 
meters  depth.  All  temperatures  were  converted  to  potential 
temperature.  Table  1  summarizes  the  remaining  data.  Each 
profile  consists  of  temperature,  salinity,  pressure,  and 
depth  measurements. 
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ITP 

Max  #  of  Obs  in 
each  Profile 

#  of  Profiles 

1 

1051 

1022 

2 

520 

122 

3 

1987 

766 

4 

1917 

349 

5 

2327 

547 

6 

1710 

667 

Table  1.  Number  of  Observations  (individual  data  points) 
for  Each  Profile.  Each  profile  has  up  to  the  maximum 
number  of  observations  indicated.  Not  all  profiles 
have  the  maximum  number  of  observations. 

These  data  were  further  reduced  by  determining  the 
values  of  temperature  and  salinity  at  the  center  of  each 

layer.  This  was  accomplished  by  searching  the  data  for  the 
local  maxima  of  vertical  temperature  gradient.  This  local 
maximum  corresponds  to  the  interface  between  each  layer. 
Using  the  position  of  each  interface,  the  center  of  each 

layer  is  determined  to  be  exactly  half  way  between  two 
adjacent  interfaces.  The  temperature,  salinity  and  depth  of 
these  data  points  is  recorded  and  marked.  Table  2  shows  the 

number  of  individual  diffusive  layers  found  in  the  data 

organized  by  profiler. 


ITP 

#  of  Layers 

#  of  Profiles 

1 

48 

1022 

2 

39 

122 

3 

51 

766 

4 

54 

349 

5 

51 

547 

6 

62 

667 

Table  2.  Number  of  Layers  for  Each  Profile. 
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Figure  6.  Temperature  -  Salinity  plot  for  ITPs  1-6. 

The  diffusive  layers  next  must  be  classified 
independently  of  the  particular  cast.  Timmermans  (2008) 
showed  that  each  diffusive  layer  can  be  identified  via  a 
plot  of  potential  temperature  versus  salinity.  Figure  6  is 
a  Temperature-Salinity  plot  for  each  ITP  used  in  this  study. 
The  vertical  clusters  at  nearly  constant  salinity  are 
associated  with  a  horizontally  homogeneous  layer,  a 
diffusive  layer.  Due  to  the  remarkable  lateral  coherence  of 
layers  in  the  thermohaline  staircase,  it  is  possible  to 
consolidate  the  data  into  one  dataset  by  identifying 
pronounced  layers  that  appear  in  all  six  datasets. 
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Figure  7.  Histogram  of  data:  Number  of  data  points  for  a 

given  salinity  value. 

Using  a  moving  window  across  salinity  space  that 
corresponds  to  a  0.001  PSU  increment  from  34.3  PSU  to  34.6 
PSU  we  estimated  the  concentration  of  data  points  as  a 
function  of  salinity,  as  shown  in  Figure  7.  The  blue  curve 
is  the  original  histogram.  The  red  curve  is  a  5-point 
moving  average,  which  is  used  to  identify  the  individual 
layers.  Noting  that  the  concentration  of  points  is  highest 
at  the  center  of  the  layers  and  lowest  at  the  interface,  we 
assume  that  the  troughs  in  the  histogram  represent  the 
boundary  (interface)  between  each  layer.  The  black  dashed 
lines  represent  the  inferred  interface  locations. 
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ITP 

#  of  Layers 

#  of  Profiles 

1 

19 

1022 

2 

19 

122 

3 

19 

766 

4 

19 

349 

5 

19 

547 

6 

19 

667 

ALL 

19 

3452 

Table  3.  Total  number  of  layers  and  profiles  in  final  data 

set . 

These  values  are  used  to  enumerate  layers  in  all 
profiles.  The  result  is  a  total  of  19  layers  common  to  each 
profiler  (Table  3)  .  The  data  between  adjacent  troughs 
(dashed  lines  in  Figure  7)  are  considered  to  be  in  the  same 
layer.  Thus,  layer  1  contains  data  only  from  the  region 
between  the  first  two  troughs  on  the  left  of  the  graph. 
Layer  1  corresponds  to  the  shallowest  layer  in  the 
thermocline,  and  layer  19  corresponds  to  the  deepest  (as 
salinity  increases  downward) .  It  is  important  to  note  that 
diffusive  layers  in  the  Beaufort  Sea  are  not  limited  to 
these  19.  These  are  the  layers  that  are  common  to  all  6 
ITPs  used  in  this  study. 

The  highlighted  portion  of  Figure  7  shows  the 
difficulty  in  capturing  each  layer  accurately.  The  blue 
line  clearly  shows  a  highly  variable  structure  in  this 
region,  and  our  method  for  determining  the  interfaces 
between  adjacent  layers  is  not  precise  enough  to  account  for 
this  variation.  As  a  result,  we  observed  some  anomalous 
results  in  the  inverse  calculations. 
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Locations  of  each  profile 


Longitude 


Figure  8.  Locations  of  each  profile. 

Ultimately,  each  profile  yielded  a  specific  value  of 
temperature,  potential  temperature,  salinity,  layer 
thickness,  depth,  pressure,  and  latitude  and  longitude  for 
each  of  the  19  layers.  Figure  8  is  a  plot  of  the  locations 
of  all  3452  profiles.  The  next  step  involves  interpolation 
of  the  data  onto  a  structured  grid  for  use  with  the  inverse 
model . 

C.  INTERPOLATING  DATA  TO  A  GRID 

The  inverse  model  used  in  this  study  requires  the  data 
to  be  applied  to  an  equidistantly  spaced  grid  in  both 
horizontal  directions.  This  is  accomplished  using  the 
MATLAB  routine  "griddata, "  which  linearly  interpolates  the 
data  to  a  pre-determined  set  of  grid  coordinates.  The 
result  is  a  dataset  that  is  compatible  with  the  inverse 
model.  The  horizontal  length  scale  is  15  km,  both  zonal  and 
meridional,  unless  otherwise  noted.  The  vertical  length 

scale  is  defined  by  the  layer  thickness  at  each  location. 
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3D  Model  Domain 


Figure  9.  Location  of  the  three-dimensional  dataset. 

The  grid  used  for  the  three  dimensional  model  has  the 
dimensions  7X7X19.  Figure  9  shows  the  location  of  the  grid. 
This  location  in  the  southern  Beaufort  Sea  was  chosen 
because  of  the  high  density  of  original  profile  locations, 
as  well  as  the  relatively  even  distribution  of  profiles 
throughout  the  grid  domain. 
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Figure  10.  (a)  Salinity  surfaces  for  19  layers  used  in  the 

inverse  calculation,  (b)  The  corresponding  potential 

Temperature  surfaces. 

The  structured  data  are  shown  in  Figure  10,  which 
reveals  nearly  constant  temperature  and  salinity  values  in 
each  of  the  diffusive  layers,  and  the  lateral  coherence  of 
those  layers  over  of  the  model  domain.  For  both  Figures  10a 
and  10b,  the  layers  are  arranged  from  the  deepest  to  the 
shallowest  from  top  to  bottom.  The  salinity  values  for  each 
layer  are  nearly  constant  spatially.  The  temperature  values 
vary  slightly,  but  the  homogeneity  is  evident. 
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III.  INVERSE  MODEL 


A.  MODEL  FRAMEWORK 


In  order  to  apply  (9),  (10)  and  (11)  to  measurements 
they  must  be  discretized  in  space.  This  is  achieved  by 
integrating  vertically  from  the  bottom  (b)  to  the  top  (t)  of 
each  layer,  assuming  vertically  uniform  temperature  and 
salinity  within  each  layer,  consistent  with  characteristic 
properties  of  diffusive  convection  in  thermohaline 
staircases . 


V* 


V* 
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\T 
V  J 
f-  \ 

ys 
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V  ) 

VDV^  +  [>v]^  =0 


(12) 

(13) 

(14) 


V  represents  the  vertically  integrated  horizontal 
velocity  along  the  layer,  and  w  represents  the  vertical 
velocity  across  the  interface.  The  vertical  influences  into 
and  out  of  each  layer  are  described  by  the  diffusive  and 
convective  terms  in  (12)  and  (13)  (Lee  et  al .  ,  1991)  . 


The  inverse  model  is  based  on  the  optimal  solution  to 
the  overdetermined  linear  system  of  equations  Ax=b. 
Literature  refers  to  A  as  the  data  matrix  and  b  as  the 
observation  vector.  We  will  follow  this  same  convention  in 
this  study.  The  solution  is  x,  a  vector  with  the  same 
dimensions  as  b  where  each  component  of  x  represents  the 
coordinates  of  b  in  terms  of  the  columns  of  A. 
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Figure  11.  The  seven-point  grid  network  for  finite 
differencing.  (After  Lee  &  Veronis,  1991) 


To  formulate  the  problem  in  terms  of  the  algebraic, 
rather  than  the  differential  system  the  governing  equations 
are  discretized  in  three  dimensions  using  centered  finite- 
differencing  on  a  seven-point  grid  (Figure  11)  .  Each  cell 
in  the  grid  has  the  measured  tracer  value,  temperature  and 
salinity,  in  the  center,  u  and  v  on  the  lateral  boundaries 
and  w.  Ft  and  Fs  on  the  vertical  boundaries. 
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The  model  discretization  is  employed  exactly  as 
described  in  Lee  and  Veronis  (1991).  Equations  (15),  (16) 
and  (17)  are  the  discretized  forms  of  the  advective 
diffusive  equations  and  the  continuity  equation 
respectively.  The  indices  i,  j,  k  represent  the  location  of 
the  tracer  for  that  grid  cell,  where  i  increases  eastward,  j 
increases  northward,  and  k  increases  downward  with  depth. 
The  ^  values  correspond  to  the  boundaries  of  the  grid  cell. 
Velocities  and  diffusive  fluxes  into  and  out  of  each  cell 
are  calculated  at  the  boundary  between  grid  cells.  Thus, 
the  coefficients  are  weighted  averages  of  thickness  and 
tracer  values  at  each  corresponding  boundary.  Slight 
variations  are  used  in  the  three-dimensional  models  as  we 
explore  different  formulations  of  flux  calculations.  These 
variations  are  described  in  Chapter  IV. 

Notice  that  the  integrated  forms  of  the  original 
equations  are  homogeneous.  Thus,  the  model  will  have  an 
exact  trivial  (zero)  solution.  To  avoid  this  unphysical 
solution  one  must  introduce  a  known,  or  estimated,  quantity 
that  can  be  used  to  populate  a  non-zero  observation  vector 
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b.  This  known  quantity  will  be  referred  to  as  the  reference 


unknown  x^.  To  illustrate  the  concept  of  a  reference 
variable  consider  a  homogeneous  system  (b=0)  of  four 
equations  with  three  unknowns. 

l+£  2-8  3+8 

4-8  2+8  1-8 
8+^  \-8  5+8 

6-8  5+^  3+8 

Let  us  assume  that  the  value  of  xi  is  known  or  can  be 
readily  estimated  with  some  degree  of  certainty.  All  of  the 
equations  can  be  divided  by  xi=Xr.  This  effectively  forces 
an  inhomogeneous  term  where  bAO  and  creates  a  non-trivial 
solution  to  the  problem. 

2-8  3+^^  f \+8 

2  +  8  1-8  (x2  I  4-8 

1-8  5  +  8  yx^l x^)  8  +  ^ 
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Lee  (1991)  showed  that  this  is  an  effective  method  to 
produce  the  inhomogeneous  terms  that  are  necessary  for  the 
inverse  calculations.  The  model,  however,  is  very  sensitive 
to  the  choice  of  reference  variable.  This  sensitivity  is 
explored  in  Chapter  IV. 

Another  essential  step  to  obtain  a  reliable  solution  to 
the  problem  is  to  ensure  that  the  system  of  equations  is 
overdetermined,  i.e.,  that  we  have  more  equations  than 
unknowns.  This  acts  to  smooth  out  the  errors  in  the  data. 
Each  grid  cell  contains  three  equations,  two  tracer 
(temperature  and  salinity)  equations,  and  continuity.  As 
written,  there  are  three  equations  and  five  unknowns,  u,  v, 

28 


w,  Fs,  and  Ft-  Two  important  assumptions  are  made  to  ensure 
that  there  are  fewer  unknowns  than  equations.  The  first 
assumption  is  that  w,  Fs  and  Ft  are  constant  at  each 
interface,  and  the  second  is  that  w  at  the  top  layer  is 
equal  to  Ekman  pumping  velocity  (We  ~  5  x  10”^  m/s)  at  the 
top  of  the  thermocline  (Yang,  2006) .  This  ensures  that  the 
system  is  overdetermined.  In  our  case,  we  have  1275 
equations  and  1073  unknowns 

Finally,  the  data  matrix  A  and  the  observation  vector  b 
are  built  using  the  coefficients  in  (15),  (16)  and  (17). 

This  is  accomplished  by  populating  a  matrix  with  the 
coefficients  of  the  model  unknowns  in  the  correct  location. 
Each  equation  is  represented  by  a  specific  row  and  each 
unknown  by  a  specific  column.  The  result  is  a  very  sparse 
matrix,  where  most  of  the  values  in  the  matrix  are  zero. 
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Equation  (20)  is  an  example  of  the  structure  of  the  3- 
dimensional  form  of  the  data  matrix  A,  and  n  corresponds  to 
the  number  of  u  unknowns.  A  similar  pattern  continues  for 
the  coefficients  of  v,  w,  Et,  and  Fs  completes  the  matrix. 
In  our  case  A  has  the  dimensions  1275  rows  (equations)  X 
1073  columns  (unknowns) . 

The  reference  unknown  and  its  corresponding 
coefficients  are  divided  through  as  shown  in  (19),  and  the 
vector  b  is  obtained  by  subtracting  the  coefficients  for  the 
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reference  variable  from  both  sides  of  the  equation.  The 
basic  framework  for  the  TLS  method  and  the  inverse 
calculation  is  now  complete. 

B.  THE  "BEST"  SOLUTION  TO  THE  OVERDETERMINED  SYSTEM 
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Figure  12.  Least  squares  and  total  least  squares  fits  of  a 
set  of  m=20  data  points  in  the  plane  and  just  one 
unknown  x.  ai  and  bi  represent  the  components  of  an 
example  data  matrix  a  and  an  observation  vector  b.  In 
the  least  squares  fit  the  error  only  exists  in  b  and  in 
the  total  least  squares  fit  errors  exist  in  both  a  and 
b.  (From  Markovsky  et  al . ,  2007) 
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Least  Squares  Fit 


Total  Least  Squares  Fit 


Figure  13.  Least  Squares  and  Total  Least  Squares  fits  of  a 
plane  A  and  vector  b.  The  plane  A  is  a  representation 
of  all  possible  values  for  Ax.  Ultimately  only  one 
vector  Ax  represents  the  closest  approximation  of  b. 
(a)  LS  minimizes  the  vector  -r,  where  -r  only 
represents  error  in  the  vector  b,  by  projecting  b  onto 
the  plane  A.  (b)  TLS  minimizes  -[E|r]  which  represents 
the  error  in  both  A  and  b,  needed  to  project  b  onto  the 
perturbed  plane  A.  Both  the  data  matrix  and  the 
observation  vector  are  corrected  to  find  an  exact 

solution . 


The  classical  method  of  Least  Squares  attempts  to  find 
a  "best"  solution  by  assuming  that  no  exact  solution  exists 
entirely  due  to  error  in  the  observations.  A  in  theory  is 
assumed  to  be  precisely  correct  and  some  small  correction  r 
exists  such  that  Ax=b+r  is  exactly  solvable.  (In  this  case 
-r  represents  the  errors  in  the  observation  vector.)  The 

solution  of  the  problem  is  a  vector  xeR"  corresponding  to 
the  minimum  correction  required  to  make  the  problem 
solvable,  i.e.  min||Ax-b||  =min||r||  .  Geometrically,  the 

solution  represents  either  minimizing  the  sum  of  the  squared 
vertical  distances  from  the  data  points  to  the  fitting  curve 
(Figure  12a),  or  the  projecting  the  data  vector  onto  the 
plane  spanned  by  the  columns  of  A  (Figure  13a) . 

The  total  least  squares  method  is  an  extension  of  the 

classical  Least  Squares  technique  that  takes  into  account 
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errors  in  both  the  data  matrix  A  and  the  observation  vector 
b  (Golub  et  al . ,  1980)  .  The  method  assumes  that  only 
(A+E)x=b+r  is  exactly  solvable.  E  represents  the  correction 
in  the  data  matrix  A  and  r  represents  the  correction  in  the 
observation  vector  b. 

Jm  n  2 

S£|[E|>-l!,i  (21) 

i=\  j=\ 

The  solution  corresponds  to  making  the  minimum 

correction  in  the  sense  that  min  IIFE  I  rill  ,  where  IHI  denotes 

b+r^R(A+E)  l|L  I  J|lf  II  IlF 

the  Frobenius  norm,  and  E  e  and  r  e  i?"  .  The  Frobenius 

norm  is  defined  as  square  root  of  the  sum  of  the  absolute 
squares  of  the  elements  of  the  matrix  (21) . 

Geometrically,  the  total  least  squares  solution  can  be 
described  as  either  minimizing  the  sum  of  the  squared 
orthogonal  distances  from  the  data  points  to  the  fitting 
curve  (Figure  12b),  or  the  projecting  the  data  vector  onto 
the  plane  spanned  by  the  columns  of  A+E  (Figure  13b) . 

The  Total  Least  Squares  and  Least  Squares  methods  are 
shown  schematically  in  Figure  12  and  Figure  13.  The 
difference  in  the  solutions  for  the  two  methods  amounts  to 
changing  both  the  slope  of  the  line  and  the  intercept 
(Figure  12),  or  adjusting  both  the  plane  A  and  vector  b 
(Figure  13)  to  find  a  solution,  represented  by  the  fit  line 
in  Figure  12  and  the  colored  vectors  on  the  plane  in  Figure 
13. 

There  are  two  other  important  characteristics  of  the 
TLS  solution  to  consider.  First,  the  method  implicitly 
assumes  that  the  values  of  the  data  matrix  A  and  observation 
vector  b  are  roughly  the  same.  In  Figure  12,  the  range  of 
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ai  and  bi  are  roughly  the  same,  approximately  -2  to  2.  This 
is  important  to  limiting  the  sensitivity  of  the  solution 
since  the  method  of  TLS  makes  absolute,  not  relative, 
changes  to  A  and  b.  If  the  range  of  each  input  matrix  is 
nearly  the  same,  then  the  norm  of  the  solution  vector  x  will 
be  approximately  1.  If  there  are  one  or  more  orders  of 
magnitude  difference  in  the  range  of  the  data  matrix  A  and 
observation  vector  b  then  it  is  likely  that  the  solution  to 
the  total  least  squares  problem  will  be  greatly  different 
from  the  least  squares  problem.  A  large  difference  between 
TLS  and  LS  can  be  interpreted  as  a  warning  sign  that  the 
solution  may  be  physically  inconsistent,  as  illustrated  in 
Figure  14.  A  difference  by  a  factor  of  four  in  the  norms  of 
A  and  b  can  result  in  dramatic  differences  in  solutions. 
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Figure  14.  The  importance  of  similar  norms  for  A  and  b.  (a) 

A  properly  scaled  problem,  where  A  and  b  have  the  same 
norm.  (b)  The  norms  of  A  and  b  differ  by  a  factor  of 

four . 

The  second  point  is  that  the  method  implicitly  assumes 
that  the  relative  errors  associated  with  the  data  matrix  A 
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and  observation  vector  b  are  comparable.  A  large  error  in 
either  input  matrix  will  skew  the  TLS  solution  relative  to 
the  pattern  realized  in  nature. 

An  accurate  and  reliable  estimation  of  E  and  r  are 
paramount  to  the  TLS  solution.  The  value  of  the  correction, 
[E|r]  ,  tells  us  how  far  to  move  or  perturb  the  space  of 
[A|b]  to  obtain  the  solution  x.  This  is  accomplished  via 
the  Singular  Value  Decomposition  (SVD) . 


This  process  begins  with  the  augmentation  of  the  data 
matrix  A  with  the  observation  vector  b  giving  a  new  matrix 
C=[A|b] .  The  matrix  C  is  factored  into  three  distinct 
components  orthonormal  matrices  U  and  V  and  a  diagonal 
matrix  E. 


C  =  [A|b]  =  U2:V'^ 
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Equations  (22),  (23),  and  (24)  are  three  formulations 

of  the  SVD.  The  result  is  a  series  (24)  where  the  first 
term  represents  the  largest  singular  value  of  the  data,  an 
analogue  to  the  largest  amount  of  energy,  and  the  last  term 
represents  the  least  amount  of  energy  in  the  system.  This 
is  very  similar  to  Fourier  series  where  the  terms  are 
multiple  harmonics  from  lower  frequencies  to  higher 
frequencies.  SVD  simply  goes  from  higher  energy  to  lower 
energy . 
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An  important  distinction  of  SVD  is  related  to  the 

ordering  of  the  coefficients  Oi .  The  magnitude  of  each 
value  of  Oi  decreases  with  i,  and  therefore  the  most  energy 
in  the  system  is  contained  in  the  first  component  and  the 

least  energy  is  contained  in  the  last  component  of  the 

system. 

As  noted  before  the  TLS  solution  corresponds  to  the 

correction  min  ||rE|rl||  where  ||•||  denotes  the  Frobenius 

b+reR<''"'^>  "  '  II  11^ 

norm.  Given  the  SVD,  the  Frobenius  norm  can  be  easily 
calculated . 

||•||/^  =  +  ^2+"'  + (25) 

The  Frobenius  norm  (25)  differs  greatly  from  the  2-norm 
used  in  the  Ordinary  Least  Squares  calculation.  Recall,  the 
solution  to  the  Least  Squares  problem,  min||Ax-b||  =min||r||  , 

where,  IjCH^  =  <Tj  .  The  2-norm  only  considers  the  singular 

value  with  the  highest  energy  (oi) ,  whereas,  the  Frobenius 
norm  considers  a  range  of  singular  values.  Thus,  the 
Frobenius  norm  contains  a  more  inclusive  representation  of 
the  system. 


C  =  +  •  •  •  + 


(n)-,(n)T 
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cr„^i  =  mm 

rank(C+A)<n+l 


(26) 

(27) 

(28) 


Given  (24)  ,  the  closest  (in  a  sense)  rank  n  C  to 
C  is  given  by  (26),  which  is  simply  (24)  with  the  last  term 
dropped.  Therefore  (27)  and  (28)  are  true. 
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X 


(29) 


Consider  C=[A|b],  which  has  rank  (n+1)  .  A  is  rank  n 
and  when  augmented  with  b  the  rank  becomes  (n+1)  .  The  TLS 
problem  can  be  reformulated  as  (29),  which  implies  that 

C  =  [A  +  E|b  +  r]  is  no  more  than  rank  n,  because  C  is  equal  to 

C  with  the  last  term  removed.  The  last  term  in  C  is  used  to 
formulate  the  correction  matrix  A. 


Ultimately,  what  is  sought  is  a  matrix  A  =  [E  |  r]  e 
that  describes  the  minimum  perturbation  to  the  data  matrix 

A,  and  the  observation  vector  b  such  that  C  is  rank  n  and 
the  solution  for  (29)  exists. 

It  follows  that  the  solution  can  be  determined  via  the 


last  (n+1)  column  of  V,  defined  as 


where  y  e  R"  , 


and  a^O.  If  x  =  -[y/a],  and  A  =  [E|r]  =  -Cw’’’  then  (29)  is 
solved  exactly. 


C.  SCALING  THE  MATRIX 


Ensuring  that  each  equation  is  treated  with  equal  value 
in  the  model  requires  that  we  normalize  each  row.  This  has 
two  purposes:  one  is  to  prevent  the  appearance  of  spurious 
solutions,  and  the  other  is  to  up  weight  the  more  reliable 
equations  (Lee,  1991) .  In  this  study,  we  divide  each  row  by 
its  Frobenius  norm  (for  vectors  the  Frobenius  norm  is 
equivalent  to  the  2-norm) ,  so  that  all  of  the  equations  have 
the  same  norm  in  the  system. 
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It  is  sensible  to  presume  that  the  continuity  equation 
is  less  prone  to  error  than  either  tracer  equation,  since 
there  are  no  tracer  values  in  continuity.  For  this  reason, 
we  multiplied  the  continuity  equation  by  10,  effectively 
stating  that  continuity  is  more  reliable  than  the  other  two. 
The  value  of  the  multiplier  is  important,  but  choosing  10  or 
100  ultimately  does  not  change  the  outcome.  We  simply  give 
continuity  a  higher  weight  (relative  importance  to  the 
problem)  than  the  tracer  equations. 

D.  PARAMETERIZATIONS  OF  THE  VERTICAL  FLUXES 

Four  models  are  discussed  in  this  work,  each  having  a 
different  formulation  of  the  fluxes  (Fs,  Ft)  . 

Fj.  =  Kj  —  (30) 
dz 

aF^=C[Rp)[a^Ty^  (31) 
aF^=B[X:[Rp){a^Ty^  (32) 

Equation  (30)  is  the  classical  form  of  flux  that  is 
dependent  on  the  gradient  of  the  tracer  (temperature  or 
salinity) .  This  form  is  explored  in  both  a  one-dimensional 
model  and  a  three-dimensional  model. 

Equations  (31)  and  (32)  are  based  on  Turner's  4/3  flux 
law  and  invokes  a  more  specific  formulation  proposed  by 
Kelley  (1990) .  Equation  (31)  is  explored  where  the 
formulation  of  C  (Rp)  is  unknown,  while  (32)  has  the  full 
form  of  C (Rp)  and  the  amplitude  is  unknown. 

In  Chapter  IV,  we  present  the  one-dimensional  model, 
which  corresponds  to  the  (30)  formulation  of  the  fluxes,  and 
the  three  three-dimensional  models:  Model  1,  Model  2,  and 
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Model  3,  which  correspond  to  the  flux  formulations  in 
equations  (30),  (31),  and  (32),  respectively. 

Recall,  the  solution  vector  xeR"  is  composed  of  n 
elements  that  represent  the  unknowns  (u,  v,  w,  Fs,  Ft)  in 
the  discretized  model.  Each  element  in  x  must  be  multiplied 
by  Xr  to  recover  the  correct  value  of  the  unknown.  The 
output  of  the  model  is  the  vector  x. 
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IV.  MODEL  RESULTS 


The  previous  chapters  described  the  framework  necessary 
to  implement  the  model.  In  this  chapter,  the  results  of  the 
four  different  models  are  presented.  As  discussed  in 
Chapter  III,  some  of  the  elements  of  the  model  framework  are 
adjusted  to  better  understand  the  physical  processes  at  work 
in  the  diffusive  staircases.  Each  model  will  be  described 
independently . 

A.  ONE-DIMENSIONAL  MODEL 

The  one-dimensional  model  was  developed  as  a  precursor 
to  the  larger,  more  dynamically  inclusive  three-dimensional 
model.  The  one-dimensional  model  allowed  us  to  understand 
how  the  model  worked  and  what  type  of  output  to  expect.  It 
also  allowed  us  to  observe  the  role  of  horizontal  advection 
in  the  flux  calculations. 
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In  this  case,  we  neglect  all  horizontal  velocity 
components  and  consider  a  system  where  vertical  advection  of 
salt  and  heat  are  balanced  by  diffusion.  Equations  (33), 
(34)  and  (35)  are  used  for  the  inverse  one-dimensional 
calculation.  Depth  is  measured  from  the  sea  surface  and 
therefore  vertical  velocity  is  positive  downward. 


Layer 

Potential 

Temperature 

(degC) 

Salinity 

(PSU) 

Layer 

Thickness 

(m) 

Layer 

Depth 

1 

-0.6533 

34.3128 

1.7605 

220.5104 

2 

-0.6065 

34.3307 

1.9772 

222.3872 

3 

-0.5589 

34.3484 

1.9700 

224.4895 

4 

-0.5080 

34.3654 

2.0280 

226.6684 

5 

-0.4513 

34.3847 

2.1539 

229.5241 

6 

-0.4074 

34.4006 

2.1856 

232.3383 

7 

-0.3687 

34.4140 

2.2638 

234.2862 

8 

-0.3484 

34.4243 

1.1092 

237.7177 

9 

-0.3084 

34.4338 

1.8727 

239.5300 

10 

-0.2530 

34.4521 

2.7462 

240.8740 

11 

-0.1792 

34.4767 

2.4911 

245.2952 

12 

-0.1191 

34.4967 

3.2346 

249.0272 

13 

-0.0747 

34.5121 

2.8019 

252.1884 

14 

-0.0334 

34.5261 

2.2499 

254.8196 

15 

0.0062 

34.5380 

2.3186 

257.1105 

16 

0.0449 

34.5506 

2.3997 

259.7826 

17 

0.0866 

34.5641 

3.5397 

262.5461 

18 

0.1358 

34.5825 

2.3701 

265.8477 

19 

0.1705 

34.5953 

2.8645 

267.9215 

Table  4.  One  dimensional  data  set  consisting  of  average 
values  of  potential  temperature,  salinity,  and  layer 
thickness  for  the  entire  basin. 


The  inputs  into  the  one-dimensional  inverse  model  are 
average  values  of  potential  temperature,  salinity  and  layer 
thickness.  The  total  set  of  3452  profiles  is  used  to  obtain 
averages.  These  data  are  shown  in  Table  4.  It  is  clear 


that  temperature  and  salinity  increase  with  depth,  while  the 

layer  thickness  ranges  from  just  over  1  meter  to  just  over  3 
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meters . 


The  average  thickness  is  approximately  2.33  meters. 


The  average  change  in  temperature  between  layers  is 
approximately  0.046°C,  and  the  average  change  in  salinity 
between  layers  is  approximately  0.16  PSU. 

Model  parameters  (w,  Kt,  Ks)  were  calculated  for  each 
interface  above  and  below  layers  2  thru  18.  Layers  1  and  19 
are  used  for  input  on  the  boundaries  only.  Thus  18  values 
of  the  unknowns  are  calculated.  Vertical  velocity  (w)  is 
found  to  be  constant,  which  is  consistent  with  the 
integrated  form  of  the  continuity  equation.  As  expected, 
heat  Flux  (Fh)  is  upward  (positive),  and  Kt  is  negative. 


Interface 

Fh(W/m''2) 

Kt(10''-5  m2s''-1) 

w(10''-7  m/s) 

Rrho 

1 

0.84 

-0.85 

-5 

6.23 

2 

0.88 

-0.98 

-5 

5.97 

3 

0.95 

-1.03 

-5 

5.29 

4 

0.84 

-1.06 

-5 

5.30 

5 

0.96 

-1.55 

-5 

5.54 

6 

1.51 

-1.93 

-5 

5.23 

7 

0.68 

-2.90 

-5 

7.57 

8 

1.18 

-1.35 

-5 

3.50 

9 

2.62 

-1.61 

-5 

4.82 

10 

0.99 

-1.50 

-5 

4.80 

Mean 

1.56 

-2.31 

-5.00 

5.04 

Table  5.  One-dimensional  Model  1  results. 


Table  5  is  the  results  from  the  first  one  dimensional 
model  run.  The  most  important  output  for  this  model  is  the 
Heat  Flux  (Fh)  .  The  mean  heat  flux  in  the  upward  direction 
is  1.56  W/m^.  This  is  a  critical  value  to  determining 
whether  diffusive  convection  contributes  to  ice  melt. 

The  second  result  is  the  calculated  value  of  the 
amplitude  of  C  (Rp)  .  This  value  is  0.0032  in  the  Kelley 
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(1990)  formulation  of  the  4/3  flux  law  (8)  .  Recall  that 
this  formulation  was  derived  using  laboratory  experiments 
that  may  not  represent  what  is  observed  in  nature.  We 
attempt  to  calibrate  this  formulation  to  best  represent  the 
Arctic  data. 
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gure  15.  Comparison  of  the  two  laboratory  representations 
of  Heat  Flux  and  the  calculated  value  of  Heat  Flux  from 
the  one-dimensional  Model  1. 


Figure  15  illustrates  that  the  heat  flux  (Fh) 
calculated  in  one-dimensional  Model  1  is  roughly  an  order  of 
magnitude  larger  for  the  Arctic  data. 


II 

(36) 

The  heat 

flux  was  calculated 

using 

(36)  . 

Where  AT 

represents  the 

temperature  change 

between 

the 

layer  above 
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and  layer  below  the  interface,  and  AZ  represents  the 
distance  between  the  centers  of  the  two  layers. 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.84 

0.13 

2 

0.88 

0.14 

3 

0.95 

0.18 

4 

0.84 

0.20 

5 

0.96 

0.14 

6 

1.51 

0.13 

7 

0.68 

0.04 

8 

1.18 

0.22 

9 

2.62 

0.22 

10 

0.99 

0.33 

Mean 

1.56 

0.18 

Table  6.  Comparison  of  heat  flux  calculated  using  the  one 
dimensional  Model  1  and  using  Kelley  (1990)  formulation 

of  the  4/3  flux  law. 


In  this  formulation  of  the  heat  flux,  the  magnitude  of 
the  flux  is  not  only  dependent  on  the  change  in  temperature 
between  the  two  layers,  but  also  on  the  size  of  the  two 
layers.  This  proposition  is  very  different  from  the  4/3 
flux  law,  which  states  that  the  heat  flux  is  independent  of 
depth  or  layer  thickness.  Comparison  with  the  earlier 
calculations  (Kelley,  1990)  reveals  an  order  of  magnitude 
difference  in  heat  fluxes  seen  in  Table  6. 
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X  10  *  Analysis  of  fit  "Kelley  (1990)  Fit"  for  dataset  "KtTLS  vs.  density  ratio" 


Figure  16.  Fit  of  C (Rp)  as  proposed  by  Kelley  (1990)  to  the 

one  dimensional  Model  1  results. 

Figure  16  represents  an  attempt  to  calibrate  the  Kelley 
(1990)  model  to  fit  the  fluxes  obtained  with  the  inverse 
model.  This  resulted  in  a  calculation  of  the  amplitude  of 
C  (Rp)  based  on  the  output  of  the  curve  fitting  routine.  The 
result  is  a  value  of  the  amplitude  of  C  (Rp)  =  0.0479,  an 
order  of  magnitude  larger  than  the  amplitude  of  C (Rp)  in  the 
Kelley  (1990)  formulation.  The  upper  and  lower  bounds,  at 
the  95%  confidence  level,  are  0.0603  and  0.0354 
respectively.  The  associated  RMSE  is  0.0119,  and  the  sum  of 
the  squared  errors  is  2.4484e-005. 

The  one-dimensional  case  suggests  that  the  4/3  flux 
laws  as  stated  by  Kelley  (1990)  and  Marmorino  and  Caldwell 
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(1976)  are  not  applicable  to  the  Arctic  diffusive  staircases 
in  their  original  form.  Instead  an  order  of  magnitude 
adjustment  must  be  made  to  the  equations  that  correspond  to 
a  new  value  of  C  (Rp)  .  This  is  similar  to  the  idea  of  a 
transfer  function  proposed  by  Wilson  (2007).  Additionally, 
the  curve  is  not  a  perfect  fit  to  the  data.  Only  five  of 
the  data  points  lie  within  the  95%  confidence  interval. 
There  is  a  much  steeper  trend  associated  with  this  data  that 
suggests  an  even  stronger  dependence  on  density  ratio.  The 
three  dimensional  model  is  more  inclusive,  and  shows  a  trend 
more  consistent  with  the  laboratory  results. 

B.  THREE-DIMENSIONAL  MODEL 

1 .  Depth  Dependent  Discretization 

Model  1  utilizes  the  full  three-dimensional 
discretization  of  the  equations  described  by  Lee  and  Veronis 
(1990) .  The  three  dimensional  data  set  described  in  Chapter 
II  is  used  as  input  to  the  model.  The  data  set  has  425 
individual  grid  cells  for  which  the  model  parameters  will  be 
calculated,  resulting  in  a  total  of  1275  equations  and  1073 
unknowns.  The  critical  values  of  parameters  (w,  Ks,  Kt) 
remain  the  same  as  the  one  dimensional  dataset  in  that  they 
are  assumed  to  be  constant  at  each  interface.  The  model 
calculated  u  and  v  components  of  velocity  as  well. 
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Interface 

Fh(W/m''2) 

Kt(10''-5  m2s''-1) 

w(10''-7  m/s) 

Rrho 

1 

2.15 

3.02 

5.00 

6.36 

2 

2.29 

2.95 

-5.88 

6.07 

3 

2.41 

2.76 

-13.58 

5.29 

4 

1.54 

2.17 

-35.91 

5.50 

5 

1.73 

2.69 

-18.59 

5.66 

6 

1.70 

3.17 

-23.94 

5.44 

7 

3.33 

2.92 

9.82 

5.52 

8 

1.36 

1.88 

-1.67 

4.34 

9 

2.27 

3.01 

29.37 

4.33 

10 

1.91 

2.88 

-6.84 

4.75 

Mean 

2.07 

2.74 

-6.22 

5.32 

Table  7.  Three  dimensional  Model  1  results. 


Table  7  is  the  results  of  three-dimensional  model  1. 
The  heat  flux  (Fh)  is  consistent  with  that  of  the  one 
dimensional  model  for  interfaces  1  through  10.  Below 
interface  10  the  heat  fluxes  begin  to  rise  to  questionable 
levels.  The  mean  heat  flux  -2.07  W/m^  is  for  the  first  10 
interfaces.  The  assumption  of  uniform  velocity  at  the 
interface  resulted  in  vastly  different  values  of  vertical 
velocity  throughout  the  data. 
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Figure  17.  Comparison  of  the  two  laboratory  representations 
of  Heat  Flux  and  the  calculated  value  of  Heat  Flux  from 
the  three-dimensional  Model  1. 


As  in  the  one  dimensional  model  it  is  important  to  show 
how  these  results  relate  to  the  4/3  flux  laws.  Figure  17  is 
a  plot  of  the  calculated  values  for  heat  flux  and  those 
using  the  laboratory  derived  formulations.  Like  the  one 
dimensional  model  calculations,  the  three-dimensional  model 
values  of  heat  flux  are  about  an  order  of  magnitude  larger 
than  those  calculated  using  the  laboratory  formulations  of 
the  4/3  flux  law. 
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Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

2.15 

0.13 

2 

2.29 

0.13 

3 

2.41 

0.15 

4 

1.54 

0.21 

5 

1.73 

0.17 

6 

1.70 

0.10 

7 

3.33 

0.07 

8 

1.36 

0.16 

9 

2.27 

0.26 

10 

1.91 

0.36 

Mean 

2.07 

0.19 

Table  8.  Comparison  of  heat  flux  calculated  using  the 
three-dimensional  Model  1  and  the  Kelley  (1990) 
formulation  of  the  4/3  flux  law. 

The  heat  fluxes  associated  with  these  values  of 
diffusivity  are  naturally  an  order  of  magnitude  larger  than 
those  calculated  using  the  Kelley  (1990)  formulation.  Table 
8  is  a  comparison  of  heat  flux  calculations.  The  heat  flux 
from  the  model  is  calculated  using  (36)  .  The  mean  value 
shown  for  the  model  calculation  is  the  mean  for  just  the 
first  10  interfaces.  The  mean  for  the  Kelley  formulation  is 
that  of  all  19  interfaces.  The  difference  between  the  two 
is  a  full  order  of  magnitude. 
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X  lo"*  Analysis  of  fit  "Kelley  (1990)"  for  dataset  "KtTLS  vs.  density  ratio" 


Figure  18.  Fit  of  C (Rp)  as  proposed  by  Kelley  (1990)  to  the 

three  dimensional  model  1  results. 

Figure  18  is  the  attempt  to  calibrate  the  4/3  flux  law. 
The  fit  curve  excludes  the  interfaces  14-18  since  the 
results  for  those  interfaces  were  extremely  far  from  the 
rest  of  the  data.  The  results  of  the  curve  fitting  were  a 
value  of  the  amplitude  of  C  (Rp)  =  0.0635.  The  upper  and 
lower  bounds  of  the  95%  confidence  interval  are  0.0805  and 
0.0464  respectively.  Nearly  all  of  the  data  points  fall 
within  these  bounds,  which  shows  that  the  general  trend  in 
the  data  is  similar  to  that  of  Kelley  (1990),  yet  the  model 
results  are  an  order  of  magnitude  larger. 
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2.  4/3  Flux  Law  Discretization  with  C(Rp)  Unknown 

In  order  to  compare  the  results  of  the  previous  section 
to  the  4/3  flux  law  more  directly,  the  model  discretization 
had  to  be  changed.  The  new  model  discretization  corresponds 
to  three-dimensional  Model  2  as  discussed  in  Chapter 
III, Section  D. 
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The  equations  remain  the  same  except  for  the  last  two 
terms.  The  coefficients  of  diffusivity  are  changed  to 
reflect  the  (AT)^^^  relationship  described  by  Turner  (1973). 
This  relationship  is  normalized  by  the  mean  layer  thickness 
and  mean  temperature  gradient  in  order  to  maintain  the  same 
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dimensionality  as  the  original  discretization.  The 
continuity  equation  remained  the  same.  The  data  set,  number 
of  equations,  and  number  of  unknowns  are  the  same  as  the 
previous  section. 

2('c(^,)(Ar)5 

Fr=  ' _ . 

2h(ATy 

2(c{R^){Asf 

F  =_i _ 

2/2(a5)^ 

The  equations  used  to  calculate  the  fluxes  are 
different  because  of  the  changes  to  the  discretization. 
Temperature  (39)  and  salinity  (40)  fluxes  are  now  a  function 
of  (AT)^^^,  and  the  fluxes  must  be  normalized  by  the  mean 
layer  thickness  and  mean  temperature  gradient  as  well. 


(39) 


(40) 


Interface 

Fh(W/m''2) 

Kt(10''-5  m2s''-1) 

w(10''-7  m/s) 

Rrho 

1 

0.69 

0.78 

5.00 

6.36 

2 

0.79 

0.98 

-12.14 

6.07 

3 

0.86 

1.04 

3.74 

5.29 

4 

0.73 

0.60 

8.01 

5.50 

5 

0.72 

0.73 

-0.90 

5.66 

6 

0.53 

0.96 

13.52 

5.44 

7 

0.45 

1.19 

-1.68 

5.52 

8 

0.56 

0.83 

3.99 

4.34 

9 

0.82 

0.74 

-8.36 

4.33 

10 

0.85 

0.50 

2.47 

4.75 

Mean 

-9.96 

-0.83 

1.37 

5.32 

Table  9.  Three-dimensional  Model  2  results. 


The  model  results  are  presented  in  Table  9.  We  again 
observe  questionable  heat  fluxes  below  layer  10,  which  we 
attribute  to  uncertainties  of  clearly  identifying  layers  11- 
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15.  Those  results  are  not  shown.  Overall  the  heat  fluxes 
are  less  than  1  W/m^.  The  mean  heat  flux  in  the  upper  10 
layers  is  0.70  W/m^. 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.69 

0.13 

2 

0.79 

0.13 

3 

0.86 

0.15 

4 

0.73 

0.21 

5 

0.72 

0.17 

6 

0.53 

0.10 

7 

0.45 

0.07 

8 

0.56 

0.16 

9 

0.82 

0.26 

10 

0.85 

0.36 

Mean 

0.70 

0.19 

Table  10.  Comparison  of  heat  flux  calculated  using  the 

three-dimensional  model  2and  the  Kelley  (1990) 
formulation  of  the  4/3  flux  law. 


Table  10  shows  that  the  fluxes 
calculation  are  a  factor  of  3-5  larger 
(1990)  calculation.  While  they  are  not 
magnitude  different  as  in  the  previous 
still  a  significant  difference. 


in  this  inverse 
than  the  Kelley 
a  full  order  of 
section,  there  is 
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C(Rrho)  vs.  Density  Ratio  (Rrho) 
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Figure  19.  Comparison  of  the  two  laboratory  representations 
of  C  (Rp)  and  the  calculated  value  of  C  (Rp)  from  the 
three-dimensional  Model  2. 


Figure  19  is  the  C (Rp)  plot.  The  general  pattern  of 
the  dif fusivities  calculated  by  the  inverse  model  is  the 
same.  However,  the  difference  between  the  model  and  the 
Kelley  (1990)  and  Marmorino  and  Caldwell  (1976)  calculations 
is  only  a  factor  of  3-5.  The  five  outliers  for  lower  values 
of  density  ratio  correspond  to  the  lower  five  layers,  which 
show  anomalous  results. 
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Analysis  of  fit  "fit  1"  for  dataset  "KtTLS  vs.  densityratio" 
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gure  20.  Fit  of  C (Rp)  as  proposed  by  Kelley  (1990)  to  the 

three-dimensional  Model  2  results. 

Removing  the  lower  five  layers  from  the  curve  fitting 
routine  allows  for  a  much  better  estimation  of  the  amplitude 
of  C  (Rp)  .  Figure  20  is  the  result  of  the  curve  fitting 
routine.  The  result  is  a  value  of  the  amplitude  of  C (Rp)  = 
0.0185.  The  upper  and  lower  bounds  are  equal  to  0.0232  and 
0.0139  respectively.  This  represents  a  significantly 
reduced  confidence  interval  compared  to  Model  1. 

3.  C(Rp)  Calculation  using  Kelley  (1990)  4/3 

Flux  Law 

Now  we  attempt  to  directly  extract  the  amplitude  of  the 

coefficient  C (Rp)  from  the  inverse  calculation.  The  entire 

Kelley  (1990)  formulation  of  the  4/3  flux  law  is  used  as  the 
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coefficient  in  the  model.  Thus,  the  result  is  an  exact 
value  for  the  amplitude  of  C  (Rp)  as  model  output.  We  assume 
that  the  pattern  of  C  (Rp)  is  correctly  captured  by  the 
laboratory  experiments,  but  the  amplitude  (B)  requires 
calibration . 
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We  also  assume  that  B  is  constant  and  that  it  is  not 
dependent  on  the  density  ratio,  so  in  the  model  it  is 
treated  as  a  single  unknown.  This  reduces  the  number  of 
unknowns  to  1038  with  the  number  of  equations  remaining 
1275.  The  last  two  terms  are  again  normalized  by  the  mean 
layer  thickness  and  mean  temperature  gradient  in  order  to 
maintain  the  same  dimensionality  as  the  original 
discretization.  The  continuity  equation  remained  the  same. 
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Interface 

Fh(W/m''2) 

w(10''-7  m/s) 

Rrho 

1 

0.49 

5.00 

6.36 

2 

0.47 

-2.53 

6.07 

3 

0.55 

-5.07 

5.29 

4 

0.79 

11.88 

5.50 

5 

0.63 

-11.01 

5.66 

6 

0.37 

23.70 

5.44 

7 

0.25 

-43.40 

5.52 

8 

0.58 

91.44 

4.34 

9 

0.96 

-108.90 

4.33 

10 

1.34 

82.35 

4.75 

Mean 

0.64 

4.35 

5.32 

Table  11.  Three-dimensional  Model  3  results. 


The  results  are  shown  in  Table  11.  The  heat  fluxes  are 
generally  below  1  W/m^,  with  the  average  heat  flux  being 

0.70  W/m^.  The  difference  is  still  significant,  and  is 

comparable  to  the  three  dimensional  model  2  results.  The 
vertical  velocities  are  also  anomalously  high  below  layer 
11 . 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.49 

0.13 

2 

0.47 

0.13 

3 

0.55 

0.15 

4 

0.79 

0.21 

5 

0.63 

0.17 

6 

0.37 

0.10 

7 

0.25 

0.07 

8 

0.58 

0.16 

9 

0.96 

0.26 

10 

1.34 

0.36 

Mean 

0.70 

0.19 

Table  12.  Comparison  of  heat  flux  calculated  using  the 

three-dimensional  Model  3  and  the  Kelley  (1990) 
formulation  of  the  4/3  flux  law. 
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Table  12  is  the  comparison  of  Model  3  results  and  the 
fluxes  calculated  using  the  Kelley  (1990)  formulation  of  the 
4/3  flux  law.  The  difference  between  the  two  is  exactly  a 
factor  of  3.717  throughout  the  data.  This  makes  sense  in 
that  all  we  changed  is  the  amplitude  of  the  calculation,  and 
thus  all  we  changed  is  the  multiplying  factor  determining 
the  fluxes. 
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Figure  21.  Comparison  of  the  two  laboratory  representations 
of  C  (Rp)  and  the  calculated  value  of  C  (Rp)  from  the 
three-dimensional  Model  3. 


Figure  21  is  the  plot  of  density  ratio  versus 
diffusivity.  The  trend  in  the  calculated  values  from  the 
inverse  model  is  the  same  as  the  laboratory  results.  The 
difference  in  magnitude  is  approximately  a  factor  of  3. 
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There  are  no  anomalous  results  since  the  whole 
parameterization  from  Kelley  (1990)  was  used. 

The  data  fitting  routine  in  this  case  is  not  necessary, 
since  the  output  of  the  inverse  calculation  is  a  direct 
calculation  of  the  coefficient  C. 
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V.  DISCUSSION  AND  CONCLUSIONS 


A.  DISCUSSION 

This  work  attempts  to  infer  the  vertical  heat  fluxes 
through  diffusive  staircases  in  the  Beaufort  Sea  using 
techniques  of  inverse  modeling.  Table  13  summarizes  the 
heat  flux  results  for  the  four  models  presented  in  Chapter 
IV.  The  values  of  the  fluxes  are  significantly  higher  than 
laboratory  extrapolation  (Kelley,  1990)  would  suggest.  The 
mean  fluxes  shown  are  for  the  upper  10  layers  in  all  of  the 
models.  The  mean  of  those  four  is  approximately  1.25  W/m^. 
This  suggests  that  the  upward  heat  flux  due  to  diffusive 
convection  is  sufficient  enough  to  contribute  to  sea  ice 
melt  in  the  southern  Beaufort  Sea. 


Interface 

Fh  (ID  Model) 

Fh  (3D  Model  1) 

Fh  (3D  Model  2) 

Fh  (3D  Model  3) 

Fh  (Kelley) 

1 

0.84 

2.15 

0.69 

0.49 

0.13 

2 

0.88 

2.29 

0.79 

0.47 

0.13 

3 

0.95 

2.41 

0.86 

0.55 

0.15 

4 

0.84 

1.54 

0.73 

0.79 

0.21 

5 

0.96 

1.73 

0.72 

0.63 

0.17 

6 

1.51 

1.70 

0.53 

0.37 

0.10 

7 

0.68 

3.33 

0.45 

0.25 

0.07 

8 

1.18 

1.36 

0.56 

0.58 

0.16 

9 

2.62 

2.27 

0.82 

0.96 

0.26 

10 

0.99 

1.91 

0.85 

1.34 

0.36 

Mean 

1.56 

2.07 

0.70 

0.64 

0.19 

Table  13.  Comparison  of  heat  fluxes  calculated  from  all 

models . 


The  most  glaring  discrepancies  in  the  calculations  are 
the  abnormally  large  flux  values  below  layer  10  in  three- 
dimensional  models  1  and  2.  One  reason  for  these  large  flux 
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values  is  the  dependence  on  layer  thickness  in  the  model 
discretization . 


Layer 

Thickness  (m) 

1 

1.78 

2 

1.59 

3 

1.54 

4 

1.80 

5 

2.12 

6 

2.24 

7 

1.90 

8 

1.14 

9 

1.56 

10 

2.38 

11 

2.45 

12 

3.22 

13 

2.27 

14 

2.67 

15 

1.99 

16 

3.36 

17 

1.72 

18 

2.91 

19 

1.53 

Table  14.  Mean  layer  thicknesses. 


Table  14  is  the  average  layer  thicknesses  for  each 
layer.  The  deeper  the  layers  the  larger  the  thickness 
values  become.  There  is  a  significant  increase  in  layer 
thickness  below  layer  10.  Layer  thickness  is  in  the 
denominator  in  the  vertical  terms.  This  increase  in 
thickness  effectively  reduces  the  value  of  the  coefficient 
in  the  discretized  model  and  requires  an  increase  in  the 
calculated  unknown.  A  discretization  scheme  that  uses  an 
average  layer  thickness  may  reduce  this  erroneous  result. 

The  one-dimensional  model  performs  well,  and  is  less 
variable  than  three-dimensional  Models  1  and  2,  and  its 
variability  is  comparable  to  three-dimensional  Model  3. 
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Neglecting  horizontal  velocity  likely  dampened  the  effect  of 
the  layer  thickness  dependence  in  the  flux  calculations. 
Thus,  there  were  no  erroneous  results  in  the  deeper  layers. 
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The  vertically  integrated  thermal  wind  equations  state 
that  horizontal  velocities  are  balanced  by  horizontal 
pressure  gradients.  This  implies  that  the  direction  and 
magnitude  of  the  horizontal  velocity  is  strictly  controlled 
by  the  thickness  differences  between  adjacent  cells  in  the 
data.  The  direction  of  the  flow  is  calculated  to  be  in  the 
direction  of  the  thinner  cell,  and  the  magnitude  is  a  result 
of  the  relative  differences  in  thickness  between  the  two 
cells . 


Layer  thickness  is  in  the  numerator  in  the  horizontal 
terms,  possibly  lowering  the  value  of  horizontal  velocity. 
It  is  likely  that  the  combination  of  the  sensitivity  of  the 
vertical  heat  flux  calculations  and  the  sensitivity  of  the 
horizontal  velocity  terms  to  the  layer  thickness  contributed 
to  the  anomalies  in  the  deeper  layers  in  the  two  models  were 
both  were  factors  in  the  calculation. 
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Model 

Amplitude  of  C(Rp) 

1D  Model 

0.0479 

3D  Model  1 

0.0635 

3D  Model  2 

0.0185 

3D  Model  3 

0.0119 

Kelley 

0.0032 

Marmorino  &  Caldwell 

0.0086 

Table  15.  Comparison  of  amplitudes  of  the  4/3  flux  law 

coefficient  C  (Rp)  . 

Table  15  compares  the  results  of  each  models  output  for 
the  coefficient  of  amplitude  of  C  (Rp)  .  As  one  may  expect, 
the  formulations  of  the  inverse  calculation  that  more 
closely  represents  the  Kelley  (1990)  formulation  of  the  4/3 
flux  law  (3D  Model  2  and  3D  Model  3)  yielded  results  that 
were  closer  to  Kelley's  value  of  the  amplitude.  However, 
these  results  are  still  a  factor  of  4-6  larger  than  Kelley's 
laboratory  results.  The  one-dimensional  model  and  the 
three-dimensional  Model  1  both  used  the  layer  thickness 
dependent  version  of  the  model  formulation  and  yielded 
results  more  than  a  full  order  of  magnitude  greater.  Caro 
(2009)  found  results  similar  to  Model  2  and  3  via  two- 
dimensional  direct  numerical  simulation,  and  Wilson  (2007) 
found  results  similar  to  the  ID  model  and  3D  Model  1 
analyzing  the  ITP  data. 

The  accuracy  of  the  calculations  is  not  directly 
addressed.  There  are  several  questions  associated  with  the 
accuracy  of  the  calculations  that  need  to  be  answered  with 
further  study.  One  very  important  measure  of  how  well  the 
TLS  algorithm  performs  is  the  condition  number.  The 
condition  number  measures  the  sensitivity  of  the  solution  of 
the  system  of  equations  to  errors  in  the  data.  A  small  (~1) 
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value  of  the  condition  number  indicates  a  well  conditioned 
system,  and  thus  a  system  that  is  not  too  prone  to  solution 
error.  The  larger  the  condition  number,  the  more  sensitive 
the  solution  is  to  errors  in  the  system.  The  condition 
number  of  the  data  matrix  A  in  the  one  dimensional  model  is 
781  indicating  that  the  system  is  quite  sensitive  to  data 
errors.  The  addition  of  more  equations  and  more  unknowns  to 
the  problem  only  increases  the  sensitivity.  The  condition 
numbers  for  three-dimensional  Models  1-3  are  on  the  order  of 
10®. 

B.  CONCLUSIONS 

The  inverse  modeling  technique  was  successful  at 
calculating  the  heat  flux  in  the  southern  Beaufort  Sea. 
There  are  sensitivities  associated  with  the  models  due  to 
the  differences  in  layer  thickness.  In  the  lower  layers, 
where  thickness  is  greater,  there  is  greater  potential  for 
large  flux  values  and  small  horizontal  velocities  due  to  the 
larger  values  of  layer  thickness.  Ultimately,  the  upper  10 
layers  showed  that  the  mean  heat  fluxes  were  1.25  W/m®,  and 
therefore  likely  to  contribute  to  sea  ice  melt. 

The  data  suggest  that  the  amplitude  of  the  exponential 
form  of  C  (Rp)  is  likely  within  the  range  of  a  factor  of  4  to 
nearly  an  order  of  magnitude  larger  than  laboratory  results 
(Kelley,  1990;  Marmorino  &  Caldwell,  1976)  indicate.  The 
application  of  the  inverse  model  has  shown  that 
extrapolation  of  laboratory  results  to  the  ocean  is  not 
perfect,  and  in  this  case,  not  representative  of  observed 
ocean  conditions.  Laboratory  experiments  provide  a 
foundation  for  which  observational  scientists  can  pose 
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hypotheses,  but  ultimately,  the  true  test  of  a  laboratory 
based  theory  is  how  well  the  observed  data  corroborates  the 
results . 

C .  RECOMMENDATIONS 

Several  questions  remain  unanswered.  First,  is  there  a 
way  to  improve  the  condition  of  the  system  of  equations? 
Second,  what  is  the  spatial  variability  of  the  heat  fluxes 
throughout  the  Beaufort  Sea?  Finally,  can  a  well  designed 
experiment,  where  data  is  collected  systematically  in  a 
grid,  improve  the  results? 

It  would  be  beneficial  to  add  a  set  of  equations  to 
constrain  velocities  to  known  physical  principles  such  as 
Thermal  Wind  or  Conservation  of  Vorticity.  Lee  and  Veronis 
(1990)  did  do  this,  and  their  results  were  more  consistent 
than  estimates  based  entirely  on  the  advection-dif fusion 
equations . 

The  spatial  variability  of  the  model  must  be  explored. 
This  study  only  used  one  location  in  the  southern  Beaufort 
Sea.  The  model  needs  to  be  run  in  several  locations  to  gain 
a  perspective  on  the  changes  in  fluxes  throughout  the 
region . 

The  possibility  of  Inverse  Modeling  using  an 
unstructured  grid  should  be  explored.  Unstructured  modeling 
is  a  rapidly  developing  field.  This  would  be  very  useful  in 
oceanography  because  of  the  nature  of  data  collection  in  the 
WHOI  ITP  program.  The  use  of  an  unstructured  grid  would 
allow  the  direct  implementation  of  these  modeling  techniques 
without  interpolation,  thus,  reducing  one  possible  source  of 
error . 
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In-situ  measurements  taken  within  several  days  would 
likely  improve  the  model.  We  recommend  that  in-situ 
measurements  need  to  be  collected  on  a  regularly  spaced  grid 
(similar  to  C-SALT  experiment  in  Tropics) .  A  clear  snapshot 
of  the  Beaufort  Sea  at  one  specific  time  would  significantly 
reduce  error  in  the  model.  The  dynamic  nature  of  the  Sea 
creates  a  source  of  error  naturally.  Modeling  of  this 
nature  utilizes  static  forms  of  the  equations  and  neglects 
any  temporal  variation.  Thus,  a  timely  set  of  in-situ 
measurements  would  likely  improve  the  results. 
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